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Polaron spectral functions are computed for heavily doped graphene-on-substrate systems using 
the diagrammatic Monte Carlo technique. The specific aim is to investigate the effects of interaction 
on ARPES spectra when the symmetry between graphene sub-lattices is broken by a substrate. 
Several polaronic features are visible, including band- flattening and changes in particle lifetimes. 
The difference between energies on each sub-lattice increases with coupling, indicating an augmented 
transport gap, while the spectral gap decreases slightly. In the absence of a gap, additional flattening 
is seen around the K point. Results demonstrate the potential difficulties of using ARPES to 
establish the size of bandgaps in the presence of substrates. 

PACS numbers: 73.22.Pr 
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I. INTRODUCTION 

The graphene hexagonal lattice leads to exceptional 
electronic properties: a zero bandgap semiconductor 
with very high mobilities, and essentially massless Dirac 
fermionW This has already led to a number of ap- 
plications, but transistors for digital applications re- 
main elusive because they require a gap. While a gap 
cannot be induced in suspended monolayer graphene, 
recent experimental work using angle resolved photo- 
emission spectroscopy (ARPES) has shown that there 
may be bandgaps in graphene on certain substrates: 
ARPES measurements have found a gap in graphene on a 
monolayer of intercalated gold on ruthenium^ and there 
has been significant debate regarding whether a gap is 
present in monolayer graphene on silicon carbid^*^. Re- 
cently, I established that such gaps may be enhanced 
by introducing strong electron-phonon coupling through 
a highly polarizable superstrate^. This paper discusses 
how ARPES measurements could be affected by pola- 
ronic effects. 

The controversy about the existence of the gap seen for 
graphene on SiC relates to the interpretation of ARPES 
measurements^. The authors of Ref. [3]claim that there is 
a gap of around 0.26eV, however it has been claimed that 
this gap is a mis-interpretation of other excitations such 
as polarons and polaritonPEl^ ]\jo work has previously 
been carried out to establish the properties of polarons 
in a gapped graphene system, and it is the purpose of this 
paper to compute spectral functions to understand how 
polaron properties manifest in ARPES measurements. 

In order to examine the band gap for graphene-on- 
substrate systems using ARPES, it is necessary to have 
states from which electrons can be ejected. This is car- 
ried out using either doping or gating to move the sys- 
tem away from half filling. For example, in Ref. 3 the 
Fermi energy, Ep is shifted by around 0.4eV. Bostwick 
et alW consider a system where Ep is shifted by 0.45eV. 
Where evidence has been found for the opening of the gap 
in the ruthenium system due to a breaking of the sym- 
metry of the two carbon sub-lattices in graphene, the 
Fermi energy was reported to be around 0.15eV below 



the Dirac point^. In this paper, I discuss the effect of 
polaron formation on the spectral functions of heavily 
gated or doped graphene. 

Effective electron-electron interactions can be induced 
via a strong interaction between the electrons in a 
graphene monolayer and phonons in a strongly polariz- 
able substrate because of limited out of plane screening, 
similar to that seen for quasi-2D materials^. Such inter- 
actions have been experimentally demonstrated between 
carbon nanotubes and a Si02 substrate^, and are neces- 
sary to account for the lower mobilities of graphene on 
Si02^^. These interactions will form polaronic states and 
affect the overall electronic structure for electrons in the 
graphene monolayer. 

In this paper, I calculate the effects of electron- 
phonon interaction on electrons in a honeycomb lattice 
(graphene) where a gap has been opened with a mod- 
ulated potential. I present results computed using the 
numerically exact diagrammatic quantum Monte Carlo 
technique, and compute spectral functions using stochas- 
tic analytic inference. Results from this paper are valid 
in the heavily doped regime, well away from half-filling. 
The paper is organized as follows: The model Hamilto- 
nian is introduced in Sec. [TTl Details of the extensions 
to DQMC specific to graphene are explained in Sec 



III 



Sec. [IV| presents detailed spectral functions for a range 
of model parameters. A summary and conclusions are 
given is Sec. |V| 



II. MODEL HAMILTONIAN 

In this section, I introduce details of the model Hamil- 
tonian used here. Three components are required to ex- 
amine polaronic interactions between graphene and sur- 
face phonons in a substrate: (a) Intersite hopping within 
the hexagonal plane, which is well known to properly ac- 
count for the band structure of monolayer graphene, (b) 
An electron-phonon interaction to account for any po- 
laronic effects from interaction with surface phonons in 
the substrate, (c) Direct Coulomb interaction between 
the electrons and substrate. 




FIG. 1. (Color online) Graphene-substrate system annotated 
with interactions and sublattices. Electron-phonon interac- 
tions between the graphene layer and substrate are poorly 
screened, and large interactions are possible. 



A complication of the Hamiltonian required to describe 
graphene, is that a basis of two atoms is needed to repre- 
sent the honeycomb lattice. This leads to a Hamiltonian 
with the form, 

Hti, = ^{(l)kaick^R.c.) (1) 

k 

kq 

^(^k-q^kih\^h_q) 

^^9kq \o^k-q^k{d\^d-q) 
kq 

+cI_gCfc(6t+5_g)] (2) 

H^^ = Y,^g{h% + d%) (3) 

i^static = ^V{ri)ni^ (4) 



[T] There are only weak interactions between electrons 
and phonons within the graphene plane, accounting for 
the very high mobility of suspended graphene. For sim- 
plicity, the out-of-plane interaction gkq is approximated 
to be momentum independent, and ^ = 0, leading to a 
local Holstein interaction. 

ilfph is the energy of the phonons in the substrate 
(treated as harmonic oscillators, and including both ki- 
netic and potential energy of the ions) , and the phonon 
dispersion ^q is also approximated to be momentum in- 
dependent. It is usual to define a dimensionless electron- 
phonon coupling, A = <l>(0, 0)/tMr^^, where the effec- 
tive interaction ^{n^n') = 2MQ^^gm[n]gm[i^']^ is 
approximated as a local interaction. Typical effects of 
the electron-phonon interaction in conventional semicon- 
ductors are the generation of polarons, leading to changes 
in the bandstructure and thus the effective mass of the 
carriers, modification of the Landau levels and changes in 
the optical properties of the material such as absorption 
peaks in the mid infra-red. 

To complete the model of graphene on a substrate, 
^static describes interaction between electrons in the 
monolayer and a static potential, V{r)^ induced by the 
substrate. Here, a modulated potential is considered 
where A sites have energy A and B sites —A leading to 
breaking of the symmetry between A and B sub-lattices 
and giving rise to a gap. Such a form has been suggested 
to explain gaps in graphene on SiC^ and graphene on 
rubidium^ . A similar electron-phonon Hamiltonian with- 
out the static potential was considered by Covaci and 
Berciu^^. 

Solution of this Hamiltonian is extremely involved, and 
complications arise because electron-phonon interactions 
are retarded. The next section describes how to solve the 
polaron problem for the graphene lattice using DQMC. 



Here, i^tb is the tight binding Hamiltonian represent- 
ing the kinetic energy of the electrons in the graphene 
monolayer, (j)k = —t'^-ex.p{ik.Si) and Si are the near- 
est neighbor vectors from A to B sub-lattices, ^i = 
a(l,V3)/2, ^2 = a(l,-V3)/2 and ^3 = (-a,0). Elec- 
trons are created on A sites with the operator aj and B 
sites with cj. 

The term i!/^ei-ph describes the electron-phonon inter- 
action. Phonons are created on A sites with 6^ and on the 
B sublattice with d]^. The Frohlich form for the electron- 
phonon interaction has been demonstrated experimen- 
tally for carbon nanotubes on Si02^, is theoretically pro- 
posed for quasi-2D systems where out of plane hopping 
is lo^^ and has also been found necessary to account for 
mobilities for graphene on substrate systems^. The in- 
teractions between electrons in the graphene and polariz- 
able ions in the substrate are shown schematically in Fig. 



III. METHOD: DIAGRAMMATIC QUANTUM 
MONTE CARLO 



I use the diagrammatic quantum Monte Carlo 
(DQMC) method to establish the properties of polarons 
on the graphene lattice^^^^. In order to take account 
of the bipartite lattice and AB modulated potential, 
the DQMC method has to be modified to include ba- 
sis. This is achieved by considering the non-interacting 
Green function to have a matrix form. A slight com- 
plication is presented by the off-diagonal terms of this 
matrix, which are complex. This could in principle lead 
to a phase problem (which is a generalized sign problem) . 
The complex phase is found to be very small when mea- 
suring the on-site Green functions of interest here. For 
imaginary time, r/ > r^ and at absolute zero, the Green 



functions are defined as follows: 



C — ( ^^^ Gab 
y Gba Gbb 

-(a(r/)at(r,)) -(a(r/)ct(r,)) 
Mrf)aKn)) -{c{Tf)c^{n)) 



(5) 



(6) 



where, 

Q^^ ^ [exp(-^B(r/ - Tj)) + exp(-£;A(r/ - n))] 

A [exp(-£;A(r/ - Tj)) - exp{-EB{rf - Tj))] 
[exp(-£;B(T/ - n)) + exp{-EA{rf - n))] 



Gbb 



Gab = 



Gba 



(8) 
A [exp(-.E;B(T/ - Tj)) - expj-EAJTf - Tj))] 



Afc [exp(-.EB(T/ - Tj)) - expj-EAJTf - n))] 



2vWF+A^ 



(9) 



i>l [eM-EBJTf - Tj)) - expj-EAJTf - n))] 



2yW\ 
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(10) 



Ea = ^1+ ^/WTW^ and Eb = ^- ^/A^T]M'■ 
Here ja allows greater control of the algorithm, but the 
true chemical potential lies at the bottom of the band. 
Thus, the results are only accurate when the electron 
density is low (i.e. the system is doped well away from 
half- filling) . For Tf < Ti and absolute zero, all Green 
functions are zero valued because the polaron only con- 
tains a single electron. Since G{Tf < Ti) = 0, vertices are 
ordered between times and r, where r is the length of 
the diagram. 

In its most basic form, the algorithm proceeds by 
inserting and removing interaction lines into or from 
the electron propagator. The propagators for the Hol- 
stein interaction with local phonons have the form 
exp(— a;oT)(5xx where X G {A^B} represents the sub- 
lattice type at the end of the propagator and 5xx is 
the Kronecker (5-function. Sub-lattice type is fixed at 
the ends of the whole diagram so that the dynamics of 
each symmetry broken sub-lattice can be probed inde- 
pendently. 

The imaginary time Green function tails off exponen- 
tially, and can vary by several orders of magnitude, which 
makes direct measurement of the Green function his- 
togram impossible within a reasonable time-scale. To 
avoid this, a Wang-Landau algorithm is used to make 
an initial guess for the histogram, so that all diagram 
lengths, r are visited a similar number of times during 
each simulation. The advantage of the Wang-Landau al- 
gorithm is that it obtains the histogram extremely fast. 
This histogram is not used directly for computation of 



the Green function, because the bin size is finite lead- 
ing to systematic errors. Rather, it is used as input for 
a reweighting procedure so that all imaginary times are 
visited^^ and the r dependent Green function is calcu- 
lated using the estimator given in Ref. .13, which corrects 
for finite histogram bin size. Proper choice of the pseudo 
chemical potential, /i, speeds up the initialization of the 
Wang-Landau algorithm. 

It is worth noting that Gab and Gba have a com- 
plex phase. A Monte Carlo procedure can be obtained 
by keeping track of this phase e*^ such that averages are 
given by {e^^w)/{e^^). For all cases considered here, the 
average phase (e*^) is found to be extremely close to 1, 
and no expectation values had a complex component af- 
ter averaging. There is no obvious reason why the phase 
should cancel (unlike in the ID case where the signs ex- 
actly cancel^^) and for large numbers of particles, the 
phase could become a problem. 

In order to obtain S2ectral functions, stochastic ana- 
lytic inference is usecP^. Green functions are built up 
from ^-functions that can be moved continuously in fre- 
quency using a separate Monte Carlo update scheme. 
The ability to construct spectral functions from contin- 
uous frequencies is necessary to obtain reliable analytic 
continuation at absolute zero where features can be very 
sharply peaked. Each configuration of the spectral func- 
tion is weighted as k; ex exp(— x^/2a) and the factor a 
is reduced from a large value until the average x^ < ^G^ 
the number of points in the Green function. Averages are 
then taken. Additional global updates (where all points 
can be shifted simultaneously) have been included in the 
procedure to ensure that the algorithm is ergodic. 



IV. RESULTS 

This section begins by examining how the opening 
of a band-gap affects the spectral functions when the 
electron-phonon coupling is switched on. The phonon en- 
ergy is set diS hQ = t throughout, which is relatively high 
so that features relating to polarons can be distinguished 
easily. Note that phonons of this energy are still in the 
adiabatic regime at around one third of the half band 
width. Naturally, this energy is much higher than that 
of any phonons in graphene, or of any surface phonons in 
the substrate. Values for A ranging from to t are also 
large for the same reason. 

Fig. |2] shows how the graphene spectral function, 
A{E)^ changes across the Brillouin zone. Spectral func- 
tions are computed from the full Green functions, Qaa{t) 
and Qbb{^)i with Q calculated on a logarithmic mesh 
with 500 points. Separate image plots can be seen for 
A and B type electrons to make the specific contribu- 
tions from each sub-lattice clear. The non-interacting 
band structure is superimposed for comparison. A mod- 
erate dimensionless electron-phonon coupling of A = 0.7 
is chosen, with the exception of panel (b), which shows 
spectral functions for the larger A = 2.8. Panels (a) and 




FIG. 2. Image plots of the graphene spectral function across the Brillouin zone. Various A and A are shown. The non- 
interacting dispersion is overlaid. At A = 0, A and B sub-lattices are symmetrical, leading to identical spectral functions, so 
only a single panel is shown for each A. Comparison with the non-interacting dispersion shows that there is some flattening of 
the band close to the K point. The spectral function is sharp close to the K point. In the vicinity of the F point, the spectral 
function is also sharp for states within an energy KQ of the bottom of the band. Quasi-particle lifetime (related to the inverse 
of the width) is greatly reduced close to the tops of the bands. Increase in A breaks the AB symmetry. Some band flattening 
is seen. There is also a weak excitation associated with B sites at higher energies, which touches the lower band. Panel (a) 
A = 0, A = 0.7, (b) A = 0, A = 2.8 (c) and (d) A = 0.5t, A = 0.7, A and B site electrons respectively and (e) and (f) A = t, 
A = 0.7, A and B site electrons. 



(b) correspond to A = 0, (c) and (d) to A = t/2 and (e) 
and (f) to A = t. 

The main features of Fig. [2] are: 

• The quasi-particle peaks are sharp at low energies 
E — Eq ^ hQ^ but broaden significantly for higher 
energies. This is especially noticeable in panels (a), 
(d) and (f). Eq is the polaron ground state energy. 

• A clearly identifiable polaron band (split off from 
the main dispersion) can be seen for large A in panel 

(b). 

• Asymmetry between electrons on site A and site B 
increases with A. 

• The spectral gap can be seen to increase with A. 
The spectral gap is slightly smaller than A due to 
broadening of the quasiparticle peak. 

• The beginnings of a flat polaron band can be seen 
for A = t in panel (f) and is just visible in panel 
(d) for A = t/2. 



At A = 0, A and B sites are symmetrical and the spec- 
tral function for each sub-lattice is identical, so only the 
spectral function for the A sub-lattice is shown in Figs. 
[2] (a) and (b). In the vicinity of the F point, the spectral 
function is sharply peaked for states within an energy of 
around hQ of the bottom of the band, and quasi-particle 
lifetime (related to the inverse of the width) is greatly 
reduced for the highest energy states close to the tops of 
the bands. This is a polaronic effect: within an energy hQ 
of the bottom of the band the electrons have insufficient 
energy to excite real phonons. At the top of the band, 
the spectral functions are broad due to interactions. 

The spectral function is also sharply peaked close to 
the K point. By examining the zero gap states (Fig. 2] 
(a) and (b)) it is possible to assess whether it is possible 
that polaronic states from this form of interaction could 
emulate a spectral gap in ARPES measurements. Exam- 
ination of Fig. [2|a) shows that there is a flattening of 
the dispersion near the K point (highlighted by compari- 
son with the non-interacting dispersion) accompanying a 
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FIG. 3. (Color online) Variation of the graphene spectral 
function across the Brillouin zone for various A and A = 0.7. 
The same data as Fig. [2] is shown with both A and B sub- 
lattices superimposed on the same plot. From top to bottom, 
A = 0, A = t/2 and A = t. The gap opening is clearly visible. 
Broadening of the spectral functions can be seen, especially 
at larger energies. 



decrease in the width of the quasi-particle peak. The flat- 
tening of the dispersion close to the K point is associated 
with a steepening of the dispersion along the KM and KF 
lines, giving the dispersion the appearance of a waterfall. 
The waterfall like features are at variance with the idea 
that polaronic features could look like a band gap'', since 
emulation of a gap would require the dispersion to show 
the opposite shape to that seen here: a very steep band 
structure close to K with a rapid change of gradient lead- 
ing to flattening between K and M. Panel (b) shows that 
the spectral weight gets broader with increased A, indi- 
cating that it is also likely that finite quasi-particle life- 
times contribute to an obscuring of any gap rather than 
any apparent gap opening. 

A major feature in panel (b) for A = 2.8 is the emer- 



gence of a fiat polaron band, which separates off from 
the lower band. Such a feature typically corresponds to 
the energy required to excite a real phonon. The large 
electron-phonon coupling leads to a significant decrease 
in the quasi-particle lifetime at large energies, however, 
the band fiattening seen around the K point (which is 
unconnected to the polaron band) persists. The ground 
state polaron energy is reduced due to the polaron self- 
interaction (seen as the offset from the non-interacting 
band at the F point). 

Increase in A breaks the symmetry between A and B 
sub-lattices, and this can be seen in Fig. ^ panels (c)- 
(f). A gap opens on increased A. A band at high en- 
ergies, roughly tracing the dispersion of the ungapped, 
non-interacting band, can be seen in panels (d) and (f), 
although its spectral weight is extremely small. The ori- 
gins of this band are unclear. Increased A also leads to 
a fiat polaron band, that separates from the main band 
(it can be seen as the lowest energy feature at the K 
point). This feature is visible in panel (d), but is partic- 
ularly clear in panel (f). This band is a consequence of 
polaron localization at large A, which increases polaron 
self-interaction (i.e. the effective A is increased by the 
localization). 

To make the gap clearer. Fig. [3] shows the same data, 
but with both spectral function types superimposed on 
the same plot. For clarity, spectral weight with a magni- 
tude less than 0.1 is not shown (variations on this order 
of magnitude are not distinguishable, and the resulting 
curve appears as a series of straight horizontal lines ob- 
scuring the plots). From top to bottom A = 0, A = t/2 
and A = t and A = 0.7 in all panels. The gap open- 
ing can clearly be seen. The effects of interaction, which 
increases the width of the spectral function at large ener- 
gies where real phonons can be created, is much clearer 
in these plots. This reduction in quasi-particle lifetime 
would be instrumental in making the gap look smaller in 
ARPES measurements. 

Figs. [4] and |5] show how the graphene spectral func- 
tion changes across the Brillouin zone when the electron- 
phonon coupling is varied and A = t (from top to bot- 
tom, A = 0.7,1.4,2.1 and 2.8). Spectral functions for 
the A sub-lattice can be seen on the left of Fig. |4] (pan- 
els (a), (c), (e), (g)) and spectral functions for the B 
sub-lattice are shown on the right hand side of the plot 
(panels (b), (d), (f) and (h)). Self-interactions lower the 
polaron energy, and it can be seen that spectral weight 
at the F point moves to lower energies relative to the 
non-interacting bands as A increases. At large A, a po- 
laron band separates from the main band, and can be 
seen as a low energy band that is almost fiat and has 
a high quasi-particle lifetime. The remnant of the non- 
interacting band is visible as a side-band above the po- 
laron band. The quasi-particle lifetime decreases dra- 
matically with increased A, seen as a broadening of the 
spectral function. Again, a shadow band can be seen 
in the spectral functions for B sites, although it is rel- 
atively weak, and decreases in weight as A is increased. 
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FIG. 4. Image plots of the graphene spectral function across Brillouin zone. A = t and various A. Panels (a) and (b) A = 0.7 
(c) and (d) A = 1.4, (e) and (f) A = 2.1 and (g) and (h) A = 2.8. A sites can be seen on the left and B on the right. A flat 
polaron band can be seen forming at an energy around KQ above the bottom of the band for very large A. The quasi-particle 
lifetime decreases dramatically with increased A. 



Examination of Fig. [s] (which shows spectral functions 
for both sub-lattices superimposed onto the same plot), 
shows that the spectral gap is robust against increase in 
A, although a number of additional excitations appear. 

Finally, Fig. [6] shows how the spectral function at the 
K point evolves as the electron-phonon coupling is in- 
creased. Three main peaks are visible: (1) The lowest 
corresponds to the polaron band, (2) electrons on the A 
sub-lattice have at least 1 excited state, and (3) B type 
electrons can be seen at the highest energies. 

In the strongly doped system, the gap between the A 
excited state and the B electron energy is essentially un- 
changed by an increase in the electron-phonon coupling, 
with the gap remaining robust. The polaron band can 
be seen splitting off from the non-interacting band, with 
the polaron energy dropping rapidly at small A, followed 
by a sustained decrease. This is related to the flattening 
of the polaron band seen in Fig. [4J 

As the electron-phonon coupling is increased, the 
quasi-particle lifetime drops, seen as a broadening of the 
peak. The increase in spectral function width indicates 
that ARPES measurements away from half-filling could 
underestimate the band gap in the presence of electron- 
phonon interactions. 

There is an increase in the energy difference between 



the bottom of the A band and the bottom of the B band 
at the Dirac point, indicating an increase in the potential 
barrier formed by the higher energy A sites. This indi- 
cates an increased transport gap at the K point from 
electron-phonon interactions consistent with perturba- 
tion theory calculations at half-filling^^ . 



V. SUMMARY AND CONCLUSIONS 

In this paper, I have used the diagrammatic quantum 
Monte Carlo technique to compute the spectral functions 
of polarons on a honeycomb lattice in the presence of 
a substrate. Extensions to DQMC were introduced to 
deal with the bi-partite graphene lattice. Spectral func- 
tions were obtained for a variety of electron-phonon cou- 
pling strengths and substrate induced sub-lattice sym- 
metry breaking using stochastic analytic inference. The 
results presented here related to heavily doped graphene. 

Electron-phonon interactions are seen to have a range 
of effects on the band structure: The quasi-particle peaks 
are sharp at low energies E — Eq < ft, but broaden signif- 
icantly for higher energies. A clearly identifiable polaron 
band forms. Asymmetry between electrons on different 
sub-lattices increases with A. 
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A gap in the spectral function is induced on increase 
of the energy difference between sub-lattices, A, but is 
slightly reduced by broadening of the quasiparticle peak 
at large A. The formation of a polaron band on increase 
of electron-phonon coupling increases the transport gap 
at the K point due to an increase in the energy differ- 
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FIG. 6. (Color online) Graphene spectral function at the 
K point. A = t and A is varied. The spectral gap (seen 
at approximately E = 5t) is robust. The spectral functions 
broaden with increased A, and the polaron band rapidly drops 
in energy leading to an enhanced transport gap. 



ence between sub-lattices. This indicates that strongly 
polarizable substrates and superstrates could be used to 
enhance transport gaps opened by a superstrate. 

Flattening of the band around the K point in the 
absence of sub-lattice symmetry breaking suggests that 
electron-phonon coupling could not emulate a gap in the 
absence of a modulated potential. Spectral gaps induced 
by substrates are seen to be reasonably robust against 
interactions at heavy doping, although a shortening of 
quasi-particle lifetime (broadening of the spectral func- 
tion) may make the gap difficult to discern. Thus, while 
the controversy regarding graphene on SiC will only be 
cleared up by further detailed measurements, the results 
calculated here indicate that care should be taken when 
using ARPES measurements to establish the absolute 
size of band gaps. 



FIG. 5. (Color online) Graphene spectral function across Bril- 
louin zone. A = t and various A. Data is as Fig. [4] but with 
results for both sub-lattices plotted together. It can be seen 
that the gap is robust against increase in A, although a num- 
ber of additional excitations appear. 
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